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We consider a mathematical model that describes the release of heparin-binding growth factors from 
an affinity-based delivery system. In the delivery system, heparin binds to a peptide which has been 
covalently cross-linked to a fibrin matrix. Growth factor in turn binds to the heparin, and growth factor 
release is governed by both binding and diffusion mechanisms, the purpose of the binding being to slow 
growth factor release. The governing mathematical model, which in its original formulation consists of 
five partial differential equations, is reduced to a system of just two equations. We identify the governing 
non-dimensional parameters that can be varied to tune the growth factor release rate. In particular, we 
identify a parameter regime that ensures slow passive release (usually desirable) of at least a fraction of 
the growth factor. It is found that slow release is assured if the matrix is prepared with the concentration of 
cross-linked peptide greatly exceeding the dissociation constant of heparin from the peptide, and with the 
concentration of heparin greatly exceeding the dissociation constant of the growth factor from heparin. 
Also, for the first time, in vitro experimental release data is directly compared with theoretical release 
profiles generated by the model. We propose that the two stage release behaviour frequently seen in 
experiments is due to an initial rapid out-diffusion of free growth factor over a diffusion time scale 
(typically days), followed by a much slower release of the bound fraction over a time scale depending on 
both diffusion and binding parameters (frequently months). 
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1. Introduction 

Background 

In verterbrates, the extracellular matrix is a comple x mixture of carbohy drates, proteins, and possi 



bly minerals, that surrounds the cells that form tissues ([Alberts et al.\ d2002l) ). The extracellular matrix 
helps cells to bind together, and regulates a number of cellular functions, such as differentiation, pro- 
liferation, migration, and adhesion. The matrix can achieve such regulation via the appropriate release 
of growth factors, for which it can act as a depot. Macromolecules within the structure of the matrix 
can bind growth factors with high affinity, enabling the matrix to serve as a growth factor reservoir. In 
response to changes in local physiological conditions (such as the occurrence of a wound, for example), 
cells may secrete enzymes that can release such growth factor depots from the matrix. This natural 
growth factor release mechanism has inspired the design of affinity-based drug delivery systems that 
mimic the retentive and protective properties of the extracellular matrix for growth factor. In this paper, 
we analyze a mathematical model that describes drug release from some such delivery systems, make 
recommendations as to how delivery system should be prepared, and, for the first time, compare the 
predictions of the model directly with experimental data. 

The extracellular matrix consists predominantly of two classes of macromolecules: glycosamino- 
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glycans, and fibrous proteins, such as collagen and fibronectin. Glycosaminoglycans are polysaccharide 
polymers that typically have a repeating unit consisting of two sugars. Heparin is glycosaminoglycan of 
the matrix that is known to bind with a number growth factors in vivo via electrostatic interactions, ex- 
amples being basic fibroblast growth factor (bFGF), nerve growth factor (NGF), neurotrophin-3 (NT-3), 
and vascular endothelial growth factor (VEGF). 

Affinity-based drug delivery systems 



Sakivama-Elbert & Hubbell (2000a) have developed a growth factor delivery system for wound 



healing that exploits heparin's ability to bind electrostatically with growth factor. The various elements 
of their system are depicted schematically in Figure [T] The natural blood clotting matrix, fibrin, was 
chosen as the base material. Three dimensional fibrin hydrogel scaffolds were fabricated, into which 
invading cells could infiltrate, and release (via enzymatic processes) growth factor attached to the fibrin 
matrix. The growth factor attaches to the matrix via a bi-domain peptide bound to heparin, as we now 
explain. The peptide contains a domain which covalently cross-links to the fibrin matrix. However, 
the peptide also contains a domain that can bind to heparin, and such bound heparin can in turn bind 
to heparin-binding growth factor. Hence, growth factor attachment to the matrix is dependent on three 
distinct interactions, which we crudely represent by: (fibrin)-(peptide)-(heparin)-(growth factor). 




i • Free growth factor (G) 

\<*&> Free heparin (H) 

: W Immobilized peptide (P) 

\t^t> G rowtn factor- 

; - heparin complex (GH) 

Heparin- peptide complex (HP) 



Growth factor- heparin- 
-peptide complex (GHP) 



Fibrin matrix 

Fig. 1. A schematic representation of the fibrin matrix containing all six species of the model. 



The peptide is susceptible to cleavage by enzymes released from invading cells. Cells infiltrate the 
scaffold, release enzymes that degrade the peptide, and thereby release growth factor. In these systems, 
it is desirable that the growth factor be retained by the matrix until such time as it is actively released by 
cells. However, there inevitably will be some passive release, whereby free growth factor diffuses out of 
the system before cells have had the oppurtunity to actively release it. Growth factor binding to heparin 
is reversible, and it will not be permanently fixed to the matrix even in the absence of cells. The passive 
rel ease of growth factor is usually unde sirable, and in this study we use a mathematical model proposed 
by Sakivama-Elbert & Hubbell (2000a) to help identify conditions that minimize it. 

Sakiyama-Elbert & Hubbe ll (2000a), the system was used to deliver the growth factor bFGF. 



In 



Specifically, they carried out experiments in which they placed dorsal root ganglia from chickens in 
fibrin matrices loaded with bFGF. The purpose of their experiments was to evaluate the effect of the 



OPTIMAL CONDITIONS FOR SLOW PASSIVE RELEASE OF GROWTH FACTORS 



3 Of ED 



delivery system on neurite extension from the dorsal root ganglia. Their results demonstrated that the 
delivery system could enhance neurite extension by up to about 100% relative to unmodified fibrin . 
In subseq uent studies, the sy stem ha s been used to deliv er NGF dSakivama-Elbert & Hubbelil feOOObl) 
and IWood et all d2007l Eoom NT -3 dTavlor et al\ d2004 and |Wjllerfh_eLa£| j2008b), glial-derived neu- 



rotrophic grow t h fact or (GDNF) (Wo od et alU 200& 



2009h 

(Will erthef all d2008l) ). and sonic hedgehog dWillerth et al\ d2008l) ). In lWood et al\ d2009). a 13 mm 



), plat elet-d erived grow t h fact or (PDGF) 



gap in a rat sciatic nerve was repaired using a silicone nerve guidance conduit containing the delivery 
system loaded with GDNF; a schematic representation of such a conduit is given in Figure [2] 



Peptide (I), Heparin (v»0 
Nerve Conduit and Growtn Factor (#) Fibrin Matrix 






Fig. 2. A schematic representation of a nerve guidance conduit containing the growth factor delivery system. 



Matrix preparation 



We identify which elements of the delivery system may be varied to tune gro wth factor release by 

briefly sketching how the fibrin matrices used in the above studies were prepared; see lSakivama-Elbert & Hubbell 
d2000al) for the technical details and references. A fibrinogen solution was prepared by mixing the 
following components in appropriate quantity to achieve desired concentrations: plasminogen-free fib- 
rinogen from pooled human plasma, calcium ions, thrombin, peptide, heparin, and growth factor. This 
polymerization mixture was placed in a well of a 24-well tissu e culture plate and was incubated un- 
der appropriate conditions dSakivama-Elbert & Hubbelil d2000al) ) for an hour. It is clear then that the 
concentrations of the peptide, heparin and growth factor in the polymerization mixture are easily var- 
ied. However, a distinction needs to be made between the peptide that cross-links to the fibrin matrix 
and the peptide that remains unbound in the gel sin ce it is the cross-linked peptide only that forms 
part of the delivery system. In Sak ivama et al. (1999) and Schense & Hubbel 1 (1999), a procedure for 
quantifying the fraction of peptide that cross-links to the matrix is described. If the unbound peptide in- 
teracts significantly with the other f ree components of the system, then account must be taken of this in 
the mathematical modelling, and in lWood et al.\ d2007l) a quite large mathematical model incorporating 
such effects is described. However, we do not incorporate free peptide in the model considered here for 
reasons we shall discuss in Section 3. 



The system is complex 



Lin & Metters (2006) comments that the mathematical model proposed bv lSakivama-Elbert & Hubbelil 



( 2000a) contains a relatively large number of parameters (we shall see that it has three diffusivities and 
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four rate constants), and that this complicates their modelling approach. Indeed, the model in its origi- 
nal formulation does contain six differential equations, three of which contain diffusion terms. A major 
goal of the current study was to simplify the governing mathematical model. We shall show that under 
conditions of typical interest, the model can be reduced to a standard system of just two coupled partial 
differential equations governing the evolution of the concentrations of the total growth factor and total 
heparin. We shall further show that under typical conditions, the release behaviour is dominated by the 
values of just two non-dimensional parameters: the ratio of the concentration of cross-linked peptide to 
the dissociation constant of heparin from peptide, and the ratio of the initial concentration of heparin 
to the dissociation constant of the growth factor from heparin. Lin further comments that theoretical 
profiles generated by the model had yet to be directly compared with experimental data. We address 
this issue in Section 4. 



2. Mathematical modelling 

2.1 Model equations 



The m athematical model which we shall now describe was first developed by Sakivama-Elbert & Hubbell 
(2000a). There are six species in all in the model and, following ISakiyama-Elbert & Hubbelll (bOOOaL 



we use the following notation: 

P is an immobile peptide covalently fixed to the fibrin matrix; 

H is a mobile free heparin molecule that can diffuse through the fibrin matrix; 

G is a mobile free growth factor molecule that can diffuse through the fibrin matrix; 

GH is a mobile growth factor-heparin complex that can diffuse through the fibrin matrix; 

HP is an immobile heparin-peptide complex fixed to the fibrin matrix; 

GHP is an immobile growth factor-heparin-peptide complex fixed to the fibrin matrix. 

In Figure[T] we schematically represent all six species in the matrix. The possible interactions between 
these various species are described by the following four chemical reactions: 

G + H^GH G + HP^GHP 

(2.1) 

H + P^HP GH + P^GHP 

Kr Kr 

where kf,k r ,Kf,KR are the rate constants as shown. The first reaction, for example, represents the 
reversible binding of a free growth factor molecule to a free heparin molecule, with association and 
dissociation rate constants kf and k r , respectively. The other three reactions are similarly interpreted; 
see Figure Q] It is noteworthy that we have assumed that the rate constants for the first and second 
reactions above are the same, and similarly for the third and fourth reactions. This implies that we are 
assuming that the association/dissociation behaviour of growth factor for heparin does not depend on 
whether the heparin is free or bound to peptide; a similar comment applies to the binding heparin to 
peptide. 

The problems that are considered in this paper are one-dimensional, an d throughout we shall denote 



the spatial variable by x and the time variable by t . Following the notation of Sa kivama-Elbert & Hubbell 
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(feOOOal) . we denote by c G (x,t) the concentration of free growth factor G at location x and time t; the no- 
tation for the concentrations of the other five species follows in an obvious fashion. In view of the 
chemical reactions (12. Il l, and our assumptions regarding the mobility of the various species, the govern- 
ing equations for the six concentrations take the form: 

dc G d 2 c a 
dc H d 2 c H 

— = L) vl — r KfC G C M + KfC G w K F C H C P -T K R C H p, 

dt dx L 

— ^GH ^ 2 ^^-/^G^H ^r^GH K F C GH C P ~T K R C GH p, V^'A^ 

(9c P 

"3T" = K F C H C P T" K R C H p K F C GH Cp + K R C GH p, 

dC ™__, , 
^ — K/^G^HP T" /^r^GHP ~T K F^H^p K R C H p, 

^ c OHP _ , 7 , _ 

3~ — ^/^G^HP ^r C GHPTKpC GH Cp K R C GHP , 

where D G , D H and Z) GH are the diffusivities for the free growth factor, free heparin, and free growth 
factor-heparin complex, respectively; the species P, HP and GHP are taken to be immobile since the 
peptide is assumed to be covalently fixed to the fibrin matrix, and consequently the equations for their 
concentrations do not contain diffusion terms. 



2.2 Boundary and initial conditions 

We choose simple boundary conditions for ( 12.21 i that allow direct comparison with available experimen- 




Fibrin matrix 



x=0 



Fig. 3. A schematic depiction of the experimental setup for the release of growth factor from the fibrin matrix. 



tal and theoretical results for growth factor release from the delivery system. We suppose that the fibrin 
matrix occupies ^ x ^ L, with x = giving the location of a container wall through which the growth 
factor cannot penetrate, and x — L denoting the interface between the matrix and an external medium 
into which growth factor releases; see Figure|3] At the container wall, we impose no-flux conditions for 
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the mobile species, so that: 



dc r , dc H <3c rl , 

-^=0, -^=0, -^ = onx = 0. (2.3) 
ox ox ox 



At the interface between the matrix and the external medium, we impose perfect sink conditions for the 
mobile species: 

c G = 0, c„ = 0, c G „ =0 on x = L. (2.4) 



The boundary conditions ( 12.31 ), ( 12.41 ) could also model in-vitro growth factor release from a nerve guide 



tube (Saki vama-Elbert & Hub bell (2000a)) occupying — L ^ x ^ L, with x = ±L giving the location of 
the ends of the tube; see Figure [2] In this context, ( 12.31 ) are interpreted as symmetry conditions on the 
centre-line of the tube. 

Equations (12.21 ) are solved subject to the following initial conditions: 

c G = c G , c H — c H , c GH — 0, c P = c p , c HP = 0, c G]IP = at t = 0, (2.5) 



where c G ,c°,c° denote the initial concentrations of growth factor, heparin, and peptide, respectively, 
in the polymerization mixture. In our modelling, we make the simplifying assumption that all of the 
peptide in the polymerization mixture crosslinks covalently to the fibrin matrix. This is not likely to 
occur in practice, but we shall show in Section 3 how our results may be modified to take account of the 
presence of free peptide in the matrix. It should also be noted that free peptide will typically clear the 
system over a period of a few days. We conclude this section by emphasising that the key results of this 
paper are quite general for the model under typical conditions, and are not strongly dependent on the 
particular choice of boundary and initial conditions made. 

The mathematical model is now complete, and consists of equations ( 12.21 ), d2.3l ), d2.4| ), ( 12.51 ). 



2.3 Model reduction 

We now show how the model may frequently be reduced to a coupled pair of partial differential equa- 
tions. 



2.3.1 Non-dimensionalisation Before giving the non-dimensionalisation, we first note that ( 12.21 ) may 
be written in the following equivalent form: 

d , , , s _ n d 2 c G d 2 c GU 

d <9 2 c„ d 2 C ali 

3~ ( C H "i" C GH C HP ~\~ C GHP ) = D H j h D G 



^ \ H I CH I "HP I *- GHP / — *^H 0*X 2 ' 

d d 2 c GH 

■V (c GH + c GHP ) = D GH -r-^- + k f c G (c H + Cm,) - k r (c Gii + c GHP ) , (2.6) 
at ox z 

3^ (^HP "i" C GHp) = K F C P (^H ~t~ Q)h) k R (^HP ~t~ ^CHp) ) 

3^ — kfC G c^p ~t- K F c GH c P (k r -\- k r )c ghp , 
Cp -\- C H p ~t~ C GH p = c p , 
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where, for example, equation ( 12. 61 k is obtained by forming H.2\ a+ ( 12.2K + ( 12. 2k . We denote the total 
concentrations of growth factor and heparin in the matrix by c G and c„, respectively, so that: 



-GHTl-GHP) l h — t-H ~T L GH ~T t-HP ~T L GHP- 



(2.7) 



Equations d2.61 h and (12.6b ? give the evolution equations for the total growth factor and heparin, and may 
be written in conservation form as: 

where 

jT-^ D —~D — i T --D—-D — (29) 
ox ox ox ox 

give the total flux of growth factor and heparin, respectively. 

We introduce non-dimensional variables as follows: 

_ X t _ C G ^-p_ ^gh 

■ l= Ti ' = 77T77; \ ' c c = ~r j c h — ~ > c p — ~ , c g 



L' (L 2 /Z) G 



rVJ f\l f-A) r>\3 

**n 



T T *T "T 

_ ^GHP _ _ ^HP -T _ -T _ ~T _ Jg ~T _ JH 



4' HP 4' G c°' H c»' JG (D OH c°/L)' 7H (D 0H c°/L)' 

to obtain the following non-dimensional form for the governing initial boundary value problem (upon 
dropping over-bars): 

d , . * d 2 c G d 2 c GH 

=T l c G + Cqh ~T CghpJ — A^g -> 9 T -v o I 
(7? <7X OX 

o , ( w , n »^n , 1 <? 2 Cgh 

n~ C H p ~r I^Cgh "r C GH pJ/ f/iiyc J — ~\ 7 



-j V^H I «-HP I V^GH I ^ GHP J I 'IWQJ — *-*H ~\ 2 ^ 

OT OX ^7h/g CfX' 



2 



d , , „ d 



Cr, 



(QjH ~t~ C g hp) — ^G -i ~t~ KfaCQ (C H "l~ Chp) (c g h ~t~ C GH p) , 

at Ox^ 

^^7. (^?H/G C HP ~T" C GHp) — KfrpCp (TJh/G^H C GH ) (fyl/G^HP H~ C GHp) ) (2. 10) 

<?c GHP 

HP C'GHp) H~ &G{K-t>pC G uCp C G hp) 7 
T?H/P | 

c P -\- ^7n/pC H p H c ghp — J- i 

T]h/g 

dc G <?c„ <?c G „ 

= 0, =0, — =0 onx = 0, 

OX <7X ox 

c G =0, c„ = 0, c GI! =0 onx=l, 

c G = l, c H = l, c GH =0, c P =l, c HP = 0, c GHP = at t = 0, 



where: 



D D r° r° 

^ G - n ' h > * /h/g - , f/iw - 

^GH ^GH t-fj l d 



o _ ^GH o _ ^GH „ _ ^/ C H „ _ %Cp 

dc -rT2' ^~772' bH ~lT> fop_ ^r 



(2.11) 
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are the governing non-dimensional parameters. The quantities K[, Ul Ki }r give a non-dimensional measure 
of the strength of retention of growth factor by the heparin, and of heparin by the peptide, respectively. 
We denote by 

k f k f 

the dissociation constants of growth factor from heparin, and of heparin from peptide, respectively, so 
that: 

KbK=-^r, K be = -^r- (2-12) 



K° 



K D 



It is noteworthy that Kf nl involves both the concentration of available binding sites for the growth factor 
and the dissociation constant of growth factor from heparin. The case ^> 1 corresponds to the 
growth factor being strongly retained by the heparin. Conversely, K/ m <C 1 corresponds to weak retention 
of growth factor by the heparin. The parameter Kf, P is similarly interpreted in the context of heparin 
retention by the peptide. 

The non-dimensional form for the total growth factor and heparin and their fluxes are given by: 



C H — ~t~ C H p ~t~ (c GH "1" Cghp) / 7]h 



Jg ~ G dx dx ' 



(2.13) 



1 dc 



H dx TJh/g dx 



2.3.2 Param eter values In Sakivama-Elb ert & Hub bell (200 0al) ; lTavlor et al.\ d2004l) : IWillerth et al. 
(2006, 120081) and lWood et all d2007l l2008h . the fibrin gels were prepared by placing 400 jx\ of poly 



merization mixture in the wells of a 24-well plate. The diameter of each well in such a plate is 1 .56 
cm, fro m which it follows that the thickness of the gels was L w 0.2 cm. In lSakivama-Elbert & Hubbeli 
(12000al) . the heparin diffusivity was taken to be D u = 3.13 x 10~ 5 cm 2 min _1 , and for bFGF, the diffu- 



sivities D G — 6.0 x 10 5 cm 2 min 1 and D G 



1.0 x 10 5 cm 2 min 1 were used. These values, which 



were based on the work of Saltzman et qZj d 19941) and lGaigalas et al.\ d 1995b . all have order of magnitude 
10~ 5 cm 2 min _1 . In lTavlor et al. I d2004 . where the growth factor being considered was NT-3, the diffu- 
sivities used were again of order 10~ 5 cm 2 min _1 . TakingD = 3.0x 10 -5 cm 2 min _1 as a representative 
diffusivity for a free species in the matrix and L — 0.2 cm, we calculate a typical diffusion time scale 
for the system to be L 2 /D ss 1 day. Hence, in a release experiment where the diffusivities are of order 
10~ 5 cm 2 min _1 , we would typically expect the unbound compone nts to clear the syst em o ver a period 



of some da ys, and this is consistent with the experimental results of iTavlor et al.\ d2004l) and lWood et al. 
d2007[l2008h . 



There are also time scales associat ed with the rate constants for the chem ica l reactions d2. 1 1>, nam ely, 
1 /k r , 1 /{kfcl), 1 / Kr , and 1 /(k f c°). In lSakivama-Elbert & Hubbelll d2000al) and lTavlor et al.\ d20o4 . the 
values of the rate constants for the binding of heparin to the peptide were taken to be k r w 80 min -1 
and Kp w 10 9 M _1 min -1 , and c° was of order 10 -4 M. For these values, we find that 1/k r w 1 s, 
1/(k f c p ) = 10 -5 s, and we note that these times are tiny compared to the typical time scales associated 
with diffusion (days), and furthermore, would remain so even if we made c° orders of magnitude smaller. 
For the binding of growth factor to heparin, the rate constants will depend on the nature of the g rowth 
factor, and data is unfortunately frequently lacking. For bFGF. ISakiyama-Elbert & Hubbell ( 2000al) use 
the values k r w 1 min -1 , kf w 10 8 M~ ! min -1 , and take c° w 6 x 10 -5 M. For these values, \/k r = 1 
min and 1 / (kfC° H ) w 10~ 2 s, and these times are also small compared to the diffusion time scales. 
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For NT-3, the kf and k r values are unknown, but iTavlor et al.\ (120041) gives the approximation 
Kp" « 10~ 6 M for the dissociation constant, which would imply that NT-3 has a low affinity for the 
heparin binding site. By contrast, bFGF has a relatively high affinity for the heparin binding site, with 
dissociation constant K^ u » 10~ 8 M. However, we shall show in this paper that, provided the govern- 
ing mathematical model is appropriate, slow passive growth factor release is achieveable even for low 
affinity binding of growth factor to heparin. 

2.3.3 Reduction to a pair of coupled partial differential equations We conclude from the remarks 
above that for many systems of practical interest, the time scales for the association and dissociation 
rates in the chemical reactions ( 12. It are much shorter than the diffusion time scales, so that we frequently 
have: 

L 2 L 2 L 2 1 1 11 

> 



D G D H D GH k r k f c 

H K R K F^p 

In such cases, diffusion is rate limiting since it is the slowest process. We restrict our attention to such 
systems in the current analysis. In terms of the dimensionless parameters d2.1 U . the conditions above 
imply that 8g <C rnin(Kj, H , 1) and 5n <C win(Kj, P , 1), and so the differential equations (12.101 h. d2.10l u and 
d2.10l )s are replaced by the algebraic expressions: 

^JH^-GV^H ~l~ £hp) = C GH ~F C GHP ) Kbv C ? [j]\i!G C H ~T~ C G\\) = ^H/G ~F Q]HP ) (2.14) 
^(^Q?H C G C HP C GHp) "F Kbp^-CH^P C GHP = 0, 

where 9 — 8a/ 8g — &,•/%. The first two equations in ( 12.141 ) correspond to the equilibrium forms for the 
binding of growth factor to heparin, and of heparin to peptide, respectively. 

We solve the six algebraic expressions ( 12. lOb ft, ( 12.131 h. ( 12.14b for the concentrations of the six 
species c G , c H , c GH , c P , c HP , c GHP in terms of the total concentration of growth factor, cj, and heparin, c\, 
to obtain expressions of the form: 

c g = c g(c G i^ h ), c p~ c p( c h)' c gh — ^gh(c g ,C h ), C ghp = C GHP (c G ,C H ), 
c av c hp(^ g i^h)' = ^hV^g'^h)' 



We do not display these expressions here, but they can be found in dA.U of the Appendix. We note 
that these formulae can be used to calculate the equilibrium concentrations for the various species prior 
to release since both cj and cj, are known at t = 0. A numerical calculation is not required to obtain 
such quantities. In Figure |4] we plot curves for the equilibrium fraction of bound growth factor prior 
to release using the formulae ( IA.lt . We also note that it is sufficient to solve for cj, and c T H for t > as 
the concentrations for c G , c H , c GH , c P , c HP , c GHP then follow immediately from ( IA.lt . This implies that we 
can replace the problem containing five differential equations given by (12.10b by the following problem 
containing just two coupled partial differential equations: 

dc T d dc T d 

— + — i T (c T c T c T c T ) = — - + — i T (c T c T c T c T ) = 

dt dx ' GX ?^Hx7 W S ' ~\ JhV^G'^HJ^G^! ^WCJ ' 

del del 

=0, — £ =° onx = 0, (2.15) 
dx dx 

c T G = 0, c T H = on x = 1 , 
cl — 1, cl = 1 at t — 0, 
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where 

f(r T r T r T r T ) — -D*— r (r T r T ) - — r (r T r T ) 
;Vr T r T r T r T 1 — _n*_ r (V T r T ) L_ A r (V T r T ) 

./HV^G'^H'^GJV^K*/ -) '"HV^G^H/ „ 3 ^-GHV^G'^H/' 

M T?H7G 

and where the expressions for c G (c T G ,cJ { ) and c GH (Cc,cj l j) are given in ( IA. lb . We note that ( 12.15b is in a 
standard form that can be readily given to a mathematical package such as MAPLE to solve. 

3. Analysis and results 

3.1 Optimal conditions for slow passive release: strongly retained heparin and growth factor 

In a medical device such as a nerve guide tube, it is frequently required to maintain growth factor in the 
device until such time as it is actively released by invading cells. In such cases, the device should be 
designed so as to minimise passive release of growth factor via diffusion. There are five dimensionless 
parameters that can in principle be independently varied in experiments to tune the system for a given 
growth factor, and these are: 

r° r° k r° r° 

n — _!1 n — lil A — _ V — H V — p 
'' H/G — r a ' Imp — o J ° - > A ta — „ a , "-bp t^ d 

<- G <- p R "-G-H %P 

We emphasise that t he parameters K%_ v a nd K R are in principle tunable since peptides with desired proper- 
ties can be designed ( I Wood et al. ( 2007 )). However, if the peptide is also fixed, only three dimensionless 
parameters can be independently varied in the experiments, one possible choice for these being 7Jh/g, T]h/p 
and K{, H . The parameters D* and D* cannot be changed in experiments as they are fixed for a given 
growth factor. The parameters 8 G and 5n are neglected here since they are typically tiny in systems of 
interest. 

In the literature to date, the emphasis has been on experime ntally varying the parameters r? H;r , an d 



Hh/p to determine o p timal conditions for slow pas sive release; seelSakivama-Elbert & Hubbelll (12000a); 



Taylor et al.\ (120041) : I Willerth et aZ.ld2008l) and lWood et al.\ d2007l 120081) . In particular, experiments have 
been carried out for very large values of the ratio T7„ /G , and quite small values for the ratio rj H/r . However, 
we now show that if one wishes to ensure slow passive release, then the key parameters to monitor are 
Kt, H and K^, rather than rj H/G and rj H/P . More precisely, we shall show that slow release of at least a 
proportion of the growth factor is assured provided K/, H ,Ki, P ;§> 1 (with the other parameters being 0(1), 
although there are other possibilities), or, equivalently: 

c°>^ G D H and c p »tf H D p . (3.1) 

We recall that Kj, H ,Kt, P 3> 1 corresponds to strong retention of both growth factor by the heparin, and of 
heparin by the peptide. Hence, if practicable, for slow release of growth factor, the matrix should usually 
be prepared with the initial concentration of heparin being much larger than the dissociation constant of 
growth factor from heparin, and the concentration of peptide covalently cross-linked to the fibrin matrix 
peptide being much larger than the dissociation constant of heparin from peptide. We now justify this 
conclusion using an asymptotic argument and by providing numerical evidence. In particular, we shall 
demonstrate numerically that growth factor release can be relatively fast if the conditions (13. Il l are not 
met even with tj h/g 3> 1 and rj vw 1. 
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3.1.1 Asymptotics: Kb H ,Kt, P 3> 1 We write K/ m = 1 /e, = jj./e and consider the limit e — >• in 
( IA. lb for fixed 0(1) values of c T c and c„, and with /I and all the remaining dimensionless parameters in 
( IA. 1 b being (9(1). The fraction of bound drug in the matrix, which we denote by / B , is given by: 



/b 



"GH I «-GHP 



*-GHP m 



some solutions for this quantity at t = are displayed in Figure|4] For clarity, we revert to dimensional 
quantities in this Section. In the limit e — > 0, we find that: 



7T if c„ < c T G and 4 < c°, 
-£ if cj, < cl and cj, > c°, 
7T if 4 > 4 and 4 > c°, 



1 



(3.2) 



if 4 > cl and cl < cl. 



There are also three narrow transition regions at the interfaces of the four regimes listed above, but we 
omit this detail since it does not contribute to the subsequent discussion. The results of ( 13.21 ) are readily 
interpreted. Take, for example, the case / B ~ 1 for cj, > c T G and c\ < c°. Since in the current limit both the 
growth factor and heparin are strongly retained, then provided there is enough heparin to accomodate 
the growth factor (c T H > c T G ) and enough peptide to accomodate the heparin (cj, < c°), all of the growth 
factor in the matrix will be bound to leading order (/„ ~ 1). This is the desired regime for slow passive 
release, as we now confirm. The other three cases are similarly interpreted. 

To gain insight into the time it would take for the growth factor to passively release from the matrix, 
we now consider the total flux of growth factor, f a . We find as e — > that: 



where: 



del 
~77 



-D. 



del . 
1 dx 



(Do 



n \ dc « 



if cl < c T G and c T H < c°, 



if cl < cl and cj, > c°, 



°™ { { ^ } # + f? '4 ) if cl > cl and c T „ > c? , 



(a(« 



'B(cl,cl) d -£ 



ifc„ >c T and c T H < c°, 



(3.3) 



A(cl,cl) 



(cl-cl) 2 



-Dan, B(c T a ,cl) 



liclcl 



(cl-cl) 2 



(cl-cl) 2 



and assuming the dimensionless form for the derivatives arising are 0(1). The point to note here is that 
the flux of growth factor for the fourth regime, cl > c T a and c„ < c°, is 0(e) smaller than that for the 
other three. 
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In view of ( 13.2b and ( 13. 3b . it is now clear that the optimal regime for slow passive release is c T H > cj 
and cj, < c°, which indicates that the polymerization mixture for the matrix should have c° a < c° < c p . 
Hence, we expand our original recommendations for matrix preparation, ( 13.1b , to the following: 

c°>tf° H , c°»tf H D p , c°<c°<c°. (3.4) 

We have a 'large' growth factor flux for the first three cases in ( 13.3b because for each of these 
regimes, there is a substantial free component that can diffuse. For example, for the first case, cjj < cj 
and cj, < c°, we have c a ~ cj — cjj. For the boundary and initial conditions of ( 12.15b , this free drug will 
clear the system to leading order on the time scale t = 0(L 2 /D as ). However, once the free component 
has cleared, the remaining bound component in the bulk will be governed by the fourth regime of (13.3b . 
and this will clear the system on the long time scale t = 0(L 2 /(eD DH )) ^> 0(L 2 /D GH ). Similar remarks 
apply to the second and third regimes in ( 13.3b . Hence, if the matrix is prepared with Kj, H ,Kt, P 3> 1, and 
if our governing model is appropriate, then we are assured that at least a fraction of the growth factor 
will release on the slow time scale t — 0(L 2 / (eD 0H )). Furthermore, we predict that almost all of the 
growth factor will release slowly if the matrix is prepared with c°/c° > 1 and c°/c° > 1; notice that it 
is not required that these fractions be large; see Figure [4] However, we should caution that if there is 
a substantial component of free peptide (which is not included in the model described here), then there 
can still be a significant amount of free growth factor that can release on a fast diffusion time scale. 

The two stage release behaviour just described has been observed in experiments (see Section 0), 
where one sometimes sees a proportion of the growth factor releasing quickly over a period of some 
days (which could correspond to free growth factor releasing on a diffusion time scale) followed by 
much slower release of the remaining fraction (which could correspond to a strongly retained bound 
component releasing on a longer time scale such as that described above). 

Incorporating free peptide in the analysis 
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We now consider the case where a substantial fraction of the peptide remains free and competes with 
the covalently bound peptide for free heparin. We assume that heparin bound to free peptide has the 
same binding behaviour for growth factor as free heparin. The essence of our results above carry over, 
as we now explain. We suppose that the conditions (13.4b hold, and that the ratio of peptide covalently 
attached to the fibrin, r, has been quantified. Then the concentration of cross-linked peptide in the 
system is rc°, and since c° > c°, the initial concentration of heparin that is bound to cross-linked peptide 
is, to leading order, rc° . Since in turn c° > c° , the initial concentration of growth factor trapped by the 
delivery system is, to leading order, rc° G . It follows that a fraction (1 — r) approximately of the growth 
factor will diffuse out of the system over a diffusion time scale. Hence, the final recommendation, which 
we add to (13.41 . is that c° G should be chosen so that rc° G is sufficiently large for the therapy to be effective. 

3.1.2 Numerical Solutions Two different procedures were used to numerically integrate the initial 
boundary value problem ( |2.15t . In one method, simple explicit time-stepping was used to update the 
values of c G , c„, with the other quantities being then updated using dA.ll ). Centred difference approxi- 
mations were used for c GXX , c KXX , c GHXX , and the no-flux conditions on x = were handled by introducing 
a fictitious line in the usual way. In the other method, the system was numerically integrated using the 
MAPLE command pdsolve/numeric , which is based on a centred implict finite difference scheme. 
Good agreement was obtained between the two schemes and with known analytical results. 

In Figure [5] we display numerical profiles for the fraction of total growth factor that has released 
from the system as a function of time over a period of a fortnight. Four of the five release profiles 
displayed in Figure [5] (a) correspond to the case £& H) £& P 3> 1, the regime we recommend for matrix 
preparation. The other parameter values are all 0(1), and can be found either on the figure or in its 
caption. The fifth profile, the top curve in the figure, corresponds to the case of no drug delivery system, 
and has been included for comparison. In this case, all of the growth factor is free in the gel, and its 
concentration is governed by the usual linear diffusion equation (see Section 13.2. U . The other four 
curves correspond to the four asymptotic regimes identified in the previous Section. We observe for 
these that the growth factor release rate becomes slow after a period of appproximately a day, and this 
is easily interpreted. The fast initial release phase corresponds to the rapid out-diffusion of free growth 
factor on a diffusion time scale; notice that this period coincides (tellingly) with the period over which 
the growth factor releases when there is no delivery system. Once the free component has substantially 
exited the system, the remaining bound fraction releases slowly over a long time scale, as previously 
explained. In all cases, the numerical results are consistent with the asymptotic predictions. 

In Figure [5] (b), we display numerical solutions to ( 12.151 ) for Kj, P , rj HIG ^S> 1, and for various values 
of Ki, H . These parameters correspond to the heparin being strongly retained by the peptide, and the 
initial concentration of heparin greatly exceeding that of the growth factor. However, we see from the 
figure that this is not sufficient to guarantee slow release of growth factor. For values K/ m = 0(1), which 
correspond to moderate retention of growth factor, the growth factor will release over a period of some 
days. This compares unfavourably with the results displayed in Figure|5](a) where both Ki„ and K^ H are 
large and we have slow release in all cases even though the values for rj aG there are only 0(1). 

3.2 Assessing the validity of the model experimentally 

The appropriateness of the proposed governing mathematical model may be assessed experimentally for 
a particular system by simply omitting components from the polymerization mixture when preparing 
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the fibrin gels. The governing mathematical model may then reduce considerably, making it easier to 
compare its predictions with experimental data and enabling parameter estimation. We suggest that such 
simpler systems should be assessed experimentally as a preliminary to consideration of the complete 
release system. We consider two such cases, and then make a brief remark concerning the full system. 



3.2.1 No heparin If in the preparation of the fibrin matrices, no heparin is added to the fibrinogen 
solution, then c„ = c GH = c HP = c GHP = 0. The dimensional form for the governing equations reduces to: 

-zr^- = on x = 0, c G = on x = L, c G = c G at t = 0, 
ox 



and this is easily solved by separating variables (Crank ( 1975)) to obtain: 

4c° D " (-1)" +1 ( -{2n-l) 2 n 2 D G t\ ((2n-\)%x 



cos 



K 2n-l F V 4L 2 / V 2L 

n=\ \ ' \ 

The total amount of growth factor released from the fibrin matrix by time t is then given by: 

L 

M(t) = Lc G - j c G (x,t)dx, 
o 

from which it follows that the fraction of the available growth factor released by time t is: 
M(t) 8 ~ 1 ( -{2n-\) 2 n 2 D Q t \ 

This expression contains only one unknown parameter, D G , and it predicts that the growth factor should 
release on the time scale t — 0(L 2 /D G ), which for the matrices described here and typical growth factors 
corresponds to a period of some days. We suggest a four day release experiment and at least three data 
points per day for the fraction of growth factor released. This release data may then be compared with 
( 13.6b . If there is a poor match between the experimental and theoretical profiles, or if one finds that 
an unreasonable value for D G must be used to obtain an acceptable fit, then the appropriateness of the 
model for the delivery system of interest is called into question. It may be that some other process not 
incorporated in the modelling is significantly affecting the release behaviour. 



3.2.2 No growth factor and peptide If both growth factor and peptide are omitted from the polymer- 
ization mixture, then c G = c P = c GH = c HP = c GHP = 0. The only surviving species in the model is free 
heparin, and its concentration is governed by an initial boundary value problem identical in structure to 
( 13.5b : simply substitute G with H in (13.5b and (13.6b . We now suggest a four day release experiment for 
the heparin with at least three data points per day for the fraction of heparin released. The experimen- 
tal and theoretical release results may then be compared, providing a second test of the validity of the 
model. If the correspondence between theory and experiment is good, then the heparin diffusivity D u 
is estimated (of course the value obtained must be consistent with previous estimates; it must have the 
correct order of magnitude). 
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3.2.3 The full system If the model passes the tests set for it in the previous two subsections, one 
may proceed to preparing the matrix with all components included. If the matrix is prepared in accor- 
dance with the recommendations ( 13.41 i. one should experimentally observe the growth factor release rate 
become slow after a period of a few days if the model is valid. 



4. Comparison with experimental data 

We now compare the theoretical release profiles generated by the model with in vitro experimental 
release data. For experimental data where there is no delivery system or no heparin, we make the com- 
parison with the model for all times using the analytical expression (13.61 >. However, for experimental 
data where the full delivery system is present, we do not attempt to compare the model results with 
experimental data in the first two days of release since free peptide may play a significant role in this 
period, and the model does not track the concentration of this species. In fact, to incorporate the effect of 
free peptide would require the inclusion of three more species in the model: free peptide, free peptide- 



heparin complex, and free peptide-heparin-growth factor complex; see lWood et al.\ d2007l) . This would 
add three reaction-diffusion equations to the governing mathematical model, and would complicate the 
analysis considerably. However, after a period of some days, all of these species should have substan- 
tially diffused out of the system since there is no mechanism to replenish them (the covalently bound 
peptide does not dissociate), and the species that then remain do form part of the model described here. 
We should say that it is not difficult to fit the model results with complete release profiles that include 
the first few days, but this would require the selection of parameters in the model that are not compat- 
ible with the experimental conditions. The numerical solutions displayed in Figure [5] (a) do have the 
qualitative character of many experimental profiles. 

In selecting parameter values for the model, we use wherever possible the values used in the exper- 
iments. We always use the same heparin to growth factor ratio rj wa and heparin to peptide ratio tj^ as 
used in the experiments. Where estimates for the diffusivities can be found, either in the experimental 
paper in question, or elsewhere, we use them. Where a value for D G is not available, we estimate it 
by fitting an experimental profile for no delivery system to the corresponding theoretical release pro- 
file (13.6b . For D au , we then select a value which has order of magnitude 10~ 5 cm 2 /min, and which is 
such that D GH < D a ,D H . The values for K[, T = c ? /K^ ? are calculated using the given values for c° and 
K° p = 8.67 x 10~ 8 M. This is probably an over-estimate since we are beginning our simulation after 
day two and unbound peptide will have been lost, reducing the value for c°. A similar remark applies to 
the values of c° and c°. However, the values for K/, P are of the order of thousands, and adjusting them 
by a factor of two or so will have very little effect on the resulting profiles. The selection of appropriate 
values for Kj,„ = c° H /K GH is a much more delicate issue though because its values range considerably 
in the experiments and the behaviour is usually strongly dependent on this value. Unfortunately, both 
numbers involved in the calculation of Kj, H are uncertain here since c° needs to be reduced as explained 
above and only order of magnitude estimates are available for K GH . Hence, we use Kb K as a fitting pa- 
rameter, but insist that it has the same order of magnitude as c° H /K GH where c° is the value given in the 
experiment and K GH is the order of magnitude estimate for this dissociation constant. 

We have chosen to compare the model results w ith experimental da ta drawn from three studies: (a) 



Taylor et all d2004 . (b) IWood et~al\ d2007l) . and (c) IWood et all d2008b : we shall subsequently refer to 



these as (a),(b),(c), and this labelling has also been used in Figure [6] where the comparison between the 
model and the experimental data is given. In (a), (b), (c) the growth factors used are NT-3, NGF and 
GDNF, respectively. In Figure [6] we also display in each case the solution (13.6b with the appropriate 
value for D G , which is the theoretical prediction when there is no delivery system or no heparin. In 
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Figure 0(a), we display experimental data with K[, H ~ 1 and K/ m ~ 5, but we do not attempt to fit this 
data other than to comment that it follows quite closely the behaviour of the solution ( 13.61 ). 

We note that the correspondence between the model and experimental data is satisfactory in all 
cases. It is clear that for Kj m 3> 1, the experimental release rates becomes slow after a few days, which 
is consistent with the theoretical prediction. For 0(1) values of K/ m (Figure [6] (a)), the experimental 
realease rates are comparable to that for no delivery system. The experimental release profiles for data 
corresponding to no delivery system or no heparin are adequately described by d3.61 >. 

5. Discussion 

We summarise our results. 

• We have shown that for conditions of typical interest, the governing mathematical model may be 
reduced to a system of just two partial differential equations, and that the release behaviour is 
frequently dominated by the values of two non-dimensional parameters. If the model is valid for 
a particular system, there will usually be slow passive release of at least a fraction of the growth 
factor if the fibrin matrices are prepared with the concentration of crosslinked peptide greatly 
exceeding the dissociation constant of heparin from peptide, and the concentration of heparin 
greatly exceeding the dissociation constant of growth factor from heparin. It is noteworthy that 
these criteria do not preclude slow release for growth factors that bind heparin with low affinity. 
We also note the value of having reliable estimates for the two dissociation constants in the system. 

• It is experimentally convenient to vary the ratios of heparin to growth factor and of heparin to 
peptide in the polymerisation mixture for the gels to determine optimal conditions for slow passive 
release. However, these ratios are not usually the key parameters, and where this strategy does 
result in slow release, we have found that it is because the binding constants have strayed into the 
regime referred to in the point immediately above. Our results indicate that the ratios of heparin 
to growth factor and of heparin to peptide in the polymerisation mixture need neither be large nor 
small for slow release. 

• For the first time, theoretical release profiles generated by the model are compared directly with 
in vitro experimental data. It is found that once the free components have cleared the system, 
the correspondence between experimental and theoretical results is satisfactory. In particular, our 
predictions concerning conditions that will give rise to slow passive release are confirmed. 

• It may be possible to partially unpick the system experimentally by simply omitting components 
in the polymerization mixture for the fibrin gels. For example, if heparin is omitted from the 
polymerization mixture, the governing mathematical model reduces to a standard linear diffusion 
equation for the growth factor, and theoretical predictions may then be readily compared with 
experimental data to help validate the model and estimate a diffusivity; see Section [3~2l 
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Appendix 

Species concentrations in terms of total growth factor and total heparin 

Solving the six algebraic expressions (12. lOb f,, ( 12.131 ) i , ( 12.141 ) for the concentrations of the six species 
c c , c H , c aH , c P , c HP , c GHP in terms of the total concentration of growth factor, c T G , and heparin, cj,, gives: 



as. j- . fr-r _ ti Um.\ 2 4- 4t 1 h/g6 g 
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c h (c , C H ) = C H — C H p (c G , C H ) — (c GH (c G , C H ) + C GHP (c G , C H ) ) / 77h/. 

Hence, in the model considered here, it is sufficient to solve for cl and cl. 
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TABLE 1 The data used to generate the theoretical curves in Figure® Unmarked data is taken from the 
paper referred to in its column. The markings on the remaining data are explained below the table. 



Symbol 


Figure 0(a) 
NT-3 fTavlor et al. (2004")") 


Figure|6](b) 
NGF (Wood et al. (20071) 


Figure |6](c) 
GDNF (Wood et al. (200811 


Units 


D G 


2.5 x 10- 5(a) 


9.7 x 10~ 5 


2.0 x 10- 5W 


cm 2 /min 


A, 


3.0 x 10- 5W 


9.1 x 10~ 5 


3.0 x 10- sW 


cm 2 /min 


A* 


1.0 x 10- 5(c) 


7.5 x 10~ 5 


1.5xl0- 5(<,) 


cm 2 /min 




1.5xl0- 7W 


0.33 x 10~ 7(a) 


1.25xlO- 7(u) 


M 


^H-P 


8.67 x 10- 8( '° 


8.67 x 10- 8W 


8.67 x 10- 8( '° 


M 




2.3 x 10~ 4 


2.5 x 10~ 4 


2.5 x 10~ 4 


M 


c" 


7.4 x 10~ 9 


7.4 x 10~ 9 


3.1 x 10~ 9 


M 



( H ) Estimated from data in lTavlor et at . ( 2()()4): IWood ~e7ai\ j2007L[20o3) . 
W jGaigalas et al\\ 19951) . 
W lsaltzman etali |l994l). 

W K R = 78min- 1 jOlson et a/.Hl98ll) ) and k b = 9.0 x 10 s M _1 min _1 iTvler-Cross et q/.Ul994[l99l) ; lKridel et q/J j 19961) ). 
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FIG. 6. Comparison of experimental and theor etical release profile s. The curves are theoretical and the s ymbols are expe rimental. 
In the figure, (a) has data for NT- 3 taken fromlTaylor et all <2004l) . (b) has data for NGF taken from I Wood et al.\ <2007l) , and (c) 
has data for GDNF taken from I Wood etal] <2008l) . The values of Ktu used to generate the theoretical curves are given on them 
and the remaining parameters used can be found in Tablefjjof the Appendix. In (a), two approximate Aj, H values for experimental 
data are also given. 



